Preparing the data

Loading the packages and data directories

options(width=100)
out=F # do we output the plots?
render=T # do we render the plots?
qc=T # do we need quality control (stair plots)?
an=T # do we run the analyses?
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(plyr)
library(reshape)
library(matrixStats)
#library(splines)
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
if(out || render){
    library(ggplot2)
    library(gridExtra)
}
if(an){
    source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))
    library(lme4)
    library(lmerTest)
    library(BayesFactor)
}

Plot variables

if(out || render){
    # theme for plotting:
    alpha <- .7
    w <- .56 # proportion width in group plots
    xLab <- 'Target Speed'
    yLab <- 'Log Contrast Threshold'
    # colLab <- expression(paste('\nTarget\nVelocity (', degree, '/s)', sep=''))
    colLab <- expression(paste('Target Eccentricity (', degree, ')', sep=''))
    # colLabType <- 'Mask Type'
    # yLim <- c(-1.15,-0.45)
    yLim <- c(-1,-0.2)
    dodge <- position_dodge(width=0.0)
}

Plot functions

if(out || render){
    themefy <- function(p) {
        p <- p + theme_bw() + 
             theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
                axis.text=element_text(size=8), axis.title=element_text(size=9),
                legend.text=element_text(size=8), legend.title=element_text(size=9),
                legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
                legend.background = element_rect(fill='transparent'),
                plot.title=element_text(face='bold'))
    }
    sumFn <- function(ss, subjStr='subj', xStr='targV', grpStr='targEcc'){
        sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
                         mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) 
        # total mean across conditions per subj
        sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) 
        sumSubj <- merge(sumSubj, sumSubjMn)
        sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
        sumSubj$seNorm <- NA
        # sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
        sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
                      mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
                      norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
        sumGrp$subj <- 'average'
        sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
        sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
        sumComb$se[is.na(sumComb$se)] <- 0
        sumComb
    }
    plotAve <- function(pss, subjStr='subj', xStr='targV', grpStr='targEcc', 
                        xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
        pss$yMin <- pss[,yStr] - pss[,seStr]
        pss$yMax <- pss[,yStr] + pss[,seStr]
        pss[,grpStr] <- factor(pss[,grpStr])
        pss[,xStr] <- factor(pss[,xStr])
        p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                                ymin='yMin', ymax='yMax')) +
            geom_point(position=dodge, size=1, alpha=alpha) + 
            geom_line(position=dodge, alpha=alpha) +
            # scale_x_continuous(breaks=c(0,.75), labels=c('0','0.75')) +
            geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
            labs(x=xlab, y=ylab, colour=collab) + ylim(yLim) +
            guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
        p <- themefy(p)
    }
    plotIndiv <- function(pss, subjStr='subj', xStr='targV', grpStr='targEcc', 
                        xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
        pss$yMin <- pss[,yStr] - pss[,seStr]
        pss$yMax <- pss[,yStr] + pss[,seStr]
        pss[,grpStr] <- factor(pss[,grpStr])
        pss[,xStr] <- factor(pss[,xStr])
        p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                                ymin='yMin', ymax='yMax')) + 
            facet_wrap( ~ subj, ncol=4) +
            geom_point(position=dodge, size=1, alpha=alpha) + 
            geom_line(position=dodge, alpha=alpha) +
            geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
            # scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
            labs(x=xlab, y=ylab, colour=collab) + ylim(c(-1.4,0)) + 
            guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
        p <- themefy(p)
    }
}

Loading the data

allDataDir <- paste(db,'Projects/mc/data_mc3/v3 - speeds',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('mc3',dataDirs)]
colsOfInt <- c('participant', 'dom', 'mcBv', 'targBv', 'targXoff2', 
               'stairStart', 'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
    print(curDir)
    curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
    curDf <- curDf[,colsOfInt]
    if(qc){
        subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
        subjStairs <- subjStairs[grep('.tsv',subjStairs)]
        for(curStairFN in subjStairs){
            curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
            curRevs <- read.table(curStair, skip=1, nrows=1)
            curIntn <- read.table(curStair, skip=4, nrows=2)
            curInfo <- readLines(curStair)
            curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1], 
                                      stairStart=curInfo[which(curInfo==" u'startContr':")+1], 
                                      maskV=curInfo[which(curInfo==" u'mcBv':")+1], 
                                      targV=curInfo[which(curInfo==" u'targBv':")+1], 
                                      targEcc=curInfo[which(curInfo==" u'targXoff2':")+1])
            curDfRevs <- cbind(curInfoCols, curRevs[,2:11])
            nTrials <- ncol(curIntn)-1
            curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
            curDfIntn$trial <- seq(1,nTrials)
            curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
            curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
            rownames(curDfIntn) <- NULL
            dfRevs <- rbind(dfRevs, curDfRevs)
            dfIntn <- rbind(dfIntn, curDfIntn)
        }
    }
    df <- rbind(df, curDf)
}
## [1] "mc3_p1_s1_2018-07-20_1041"
## [1] "mc3_p10_s1_2018-07-23_1646"
## [1] "mc3_p2_s1_2018-07-20_1245"
## [1] "mc3_p3_s1_2018-07-20_1335"
## [1] "mc3_p4_s1_2018-07-20_1445"
## [1] "mc3_p5_s1_2018-07-20_1554"
## [1] "mc3_p6_s1_2018-07-20_1647"
## [1] "mc3_p7_s1_2018-07-23_1052"
## [1] "mc3_p8_s1_2018-07-23_1247"
## [1] "mc3_p9_s1_2018-07-23_1535"
if(qc){
    dfIntn$maskV <- round(as.numeric(levels(dfIntn$maskV))[dfIntn$maskV] * 60 / 35, 1)
    dfIntn$maskV[dfIntn$maskV<0.05] <- 0
    head(dfIntn)
}
##   subj dom stairStart maskV targV targEcc trial intn resp
## 1    1   1         -2     0  0.01      32     1 -2.0    0
## 2    1   1         -2     0  0.01      32     2 -1.5    0
## 3    1   1         -2     0  0.01      32     3 -1.0    0
## 4    1   1         -2     0  0.01      32     4 -0.5    0
## 5    1   1         -2     0  0.01      32     5  0.0    1
## 6    1   1         -2     0  0.01      32     6 -0.5    1
ds <- rename(df, c(participant='subj', meanRev6='thresh', mcBv='maskV', targXoff2='targEcc',
                   targBv='targV'))
## Compiling a data set for reversals (for analyses)
df_revs <- melt(dfRevs, id.vars=c(1:7)) # no need to specify col name strings
df_revs$variable <- as.numeric(df_revs$variable) # converting reversal names to numbers
df_revs <- rename(df_revs, c(variable='rev_n', value='rev_val'))
head(df_revs)
##   subj dom stairStart maskV targV targEcc   V2 rev_n rev_val
## 1    1   1         -2  0.01  0.01      32  0.0     1    -1.0
## 2    1   1         -2  0.01   9.6      32 -0.5     1    -1.5
## 3    1   1         -2  0.01   0.6      32 -0.5     1    -1.0
## 4    1   1         -2  0.01  0.01      64 -0.5     1    -1.0
## 5    1   1         -2  0.01   0.6      64  0.0     1    -1.0
## 6    1   1         -2  0.01   9.6      64 -0.5     1    -1.0
## Adding subject dominance thresholds:
dom_ds <- read.xlsx(paste0(db,'Projects/mc/mc3/subjs-mc3.xlsx'), 2, colIndex=c(1:8),
                    rowIndex=c(1:(length(unique(ds$subj))+1))) # 2nd sheet
ds <- merge(ds, dom_ds[,c('threshL','threshR')])
## Fixing the variable scales & factors:
ds$targEcc <- round(ds$targEcc / 35,1)
ds$maskType <- ''
ds$maskType[ds$maskV==0.01] <- 'stationary'
ds$maskType[ds$maskV==0.6] <- 'slow'
ds$maskType[ds$maskV==9.6] <- 'fast'
ds$maskV <- round(ds$maskV * 60 / 35, 1)
ds$maskV[ds$maskV<0.05] <- 0
head(ds)
##   subj dom maskV targV targEcc stairStart     thresh threshL threshR   maskType
## 1    1   1   1.0  9.60     1.8          0 -0.6783333    -1.2 -1.2005       slow
## 2    1   1   1.0  0.60     2.7         -2 -0.4216667    -1.2 -1.2005       slow
## 3    1   1   0.0  9.60     2.7          0 -0.6516667    -1.2 -1.2005 stationary
## 4    1   1   0.0  9.60     2.7         -2 -0.7216667    -1.2 -1.2005 stationary
## 5    1   1  16.5  9.60     2.7          0 -0.4183333    -1.2 -1.2005       fast
## 6    1   1   0.0  0.01     1.8          0 -0.5783333    -1.2 -1.2005 stationary

Transformed data sets

Collapsing across low and high stair starts to derive the averages for each condition.

thresh <- ddply(ds, .(subj,dom,maskV,maskType,targV,targEcc), summarise, 
                threshMean = mean(thresh))
head(thresh)
##   subj dom maskV   maskType targV targEcc threshMean
## 1    1   1     0 stationary  0.01     0.9 -0.5758333
## 2    1   1     0 stationary  0.01     1.8 -0.6333333
## 3    1   1     0 stationary  0.01     2.7 -0.4033333
## 4    1   1     0 stationary  0.60     0.9 -0.5966667
## 5    1   1     0 stationary  0.60     1.8 -0.6608333
## 6    1   1     0 stationary  0.60     2.7 -0.5900000

QC

For quality control, it’s not necessary to include the description of the conditions in the plots, as the convergence of the staircases should be the main focus.

if(qc){
    plotQc <- function(pss){
        plot(ggplot(pss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
            facet_wrap(subj ~ targV, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
    }
    a <- dlply(dfIntn, .(maskV), plotQc) 
}

Individual plots

if(render){
    plotIndivFn <- function(pss){
        ss <- sumFn(pss)
        plotSs <- ss[ss$subj!='average',]
        plot(plotIndiv(plotSs) + ggtitle(paste('maskV: ', unique(pss$maskType))))
    }
    a <- dlply(thresh, .(maskType), plotIndivFn)
}

Group plots

if(render || out){
    groupPlot <- function(pss, maskVstr){
        ss <- sumFn(pss[(pss$maskType==maskVstr),])
        themefy(plotAve(ss[ss$subj=='average',]))
    }
    p <- grid.arrange(groupPlot(thresh,'stationary') + ggtitle('stationary mask') +
                          theme(legend.position='none'),
                      groupPlot(thresh,'slow') + ggtitle('slow mask') + 
                          theme(legend.position='none', axis.title.y=element_blank()),
                      groupPlot(thresh,'fast') + ggtitle('fast mask') +
                          theme(axis.title.y=element_blank()),
                      ncol=3, widths=c(1.1,1,1.75))
    # if(render){plot(p)}
    if(out){
        jpeg('mc3.jpg', width=7, height=2.6, units='in', res=600)
        plot(p)
        dev.off()
        }
}

if(render || out){
    groupPlot <- function(pss, maskVstr){
        ss <- sumFn(pss[(pss$maskType==maskVstr),])
        themefy(plotAve(ss[ss$subj=='average',], xStr='targEcc', grpStr='targV',
                        xlab=colLab, collab=xLab))
    }
    p <- grid.arrange(groupPlot(thresh,'stationary') + ggtitle('stationary mask') +
                          theme(legend.position='none'),
                      groupPlot(thresh,'slow') + ggtitle('slow mask') + 
                          theme(legend.position='none', axis.title.y=element_blank()),
                      groupPlot(thresh,'fast') + ggtitle('fast mask') +
                          theme(axis.title.y=element_blank()),
                      ncol=3, widths=c(1.1,1,1.75))
    # if(render){plot(p)}
    if(out){
        jpeg('mc3-2.jpg', width=7, height=2.6, units='in', res=600)
        plot(p)
        dev.off()
        }
}

Analyses

## Scaling function
cent <- function(v){
    v <- apply(v,2,function(x){
        x <- x - mean(unique(x),na.rm=T)
        x <- x / max(x)
    })
}
## Centered data set
dsc <- ds
centCols <- c('dom','stairStart','targV','targEcc','maskV')
dsc[,centCols] <- cent(dsc[,centCols])
dsc$targEcc_ <- 'central'
dsc$targEcc_[dsc$targEcc==0] <- 'near'
dsc$targEcc_[dsc$targEcc==1] <- 'far'
dsc$targEcc_ <- factor(dsc$targEcc_)
head(dsc)
##   subj dom     maskV      targV targEcc stairStart     thresh threshL threshR   maskType targEcc_
## 1    1   1 -0.453125  1.0000000       0          1 -0.6783333    -1.2 -1.2005       slow     near
## 2    1   1 -0.453125 -0.4523938       1         -1 -0.4216667    -1.2 -1.2005       slow      far
## 3    1   1 -0.546875  1.0000000       1          1 -0.6516667    -1.2 -1.2005 stationary      far
## 4    1   1 -0.546875  1.0000000       1         -1 -0.7216667    -1.2 -1.2005 stationary      far
## 5    1   1  1.000000  1.0000000       1          1 -0.4183333    -1.2 -1.2005       fast      far
## 6    1   1 -0.546875 -0.5476062       0          1 -0.5783333    -1.2 -1.2005 stationary     near
## Binarized data set
dsb <- dsc
dsb$targV <- (dsb$targV + 1) / 2
dsb$maskV <- (dsb$maskV + 1) / 2
dsb$targEcc <- (dsb$targEcc + 1) / 2
head(dsb)
##   subj dom     maskV     targV targEcc stairStart     thresh threshL threshR   maskType targEcc_
## 1    1   1 0.2734375 1.0000000     0.5          1 -0.6783333    -1.2 -1.2005       slow     near
## 2    1   1 0.2734375 0.2738031     1.0         -1 -0.4216667    -1.2 -1.2005       slow      far
## 3    1   1 0.2265625 1.0000000     1.0          1 -0.6516667    -1.2 -1.2005 stationary      far
## 4    1   1 0.2265625 1.0000000     1.0         -1 -0.7216667    -1.2 -1.2005 stationary      far
## 5    1   1 1.0000000 1.0000000     1.0          1 -0.4183333    -1.2 -1.2005       fast      far
## 6    1   1 0.2265625 0.2261969     0.5          1 -0.5783333    -1.2 -1.2005 stationary     near
# pvalfn(lmer(thresh ~ dom + stairStart + maskV * targV * targEcc +
# (1|subj), data=dsb))
# pvalfn(lmer(thresh ~ dom * stairStart * maskV * targV * targEcc +
# (1|subj), data=dsb))
# pvalfn(lmer(thresh ~ dom + stairStart + maskV * targV * targEcc_b * maskSz +
#                  (1|subj), data=dsNear))

On reversals

## Centered data set
rev_c <- df_revs[df_revs$rev_n>4,]
rev_c$stairStart <- as.numeric(rev_c$stairStart)
rev_c$maskV <- as.numeric(rev_c$maskV)
rev_c$targV <- as.numeric(rev_c$targV)
rev_c$targEcc <- as.numeric(rev_c$targEcc)
centCols <- c('dom','stairStart','targV','targEcc','maskV')
rev_c[,centCols] <- cent(rev_c[,centCols])
rev_c$targEcc_ <- 'central'
rev_c$targEcc_[rev_c$targEcc==0] <- 'near'
rev_c$targEcc_[rev_c$targEcc==1] <- 'far'
rev_c$targEcc_ <- factor(rev_c$targEcc_)
rev_c$targEcc_ <- factor(rev_c$targEcc_, c('central','near','far'))
rev_c$targV_ <- 'stat'
rev_c$targV_[rev_c$targV==0] <- 'slow'
rev_c$targV_[rev_c$targV==1] <- 'fast'
rev_c$targV_ <- factor(rev_c$targV_, c('stat','slow','fast'))
rev_c$maskV_ <- 'stat'
rev_c$maskV_[rev_c$maskV==0] <- 'slow'
rev_c$maskV_[rev_c$maskV==1] <- 'fast'
rev_c$maskV_ <- factor(rev_c$maskV_, c('stat','slow','fast'))
head(rev_c)
##      subj dom stairStart maskV targV targEcc   V2 rev_n rev_val targEcc_ targV_ maskV_
## 2161    1   1         -1    -1    -1      -1  0.0     5    -0.9  central   stat   stat
## 2162    1   1         -1    -1     0      -1 -0.5     5    -0.9  central   slow   stat
## 2163    1   1         -1    -1     1      -1 -0.5     5    -0.7  central   fast   stat
## 2164    1   1         -1    -1    -1       0 -0.5     5    -0.9     near   stat   stat
## 2165    1   1         -1    -1     1       0  0.0     5    -1.0     near   fast   stat
## 2166    1   1         -1    -1     0       0 -0.5     5    -1.0     near   slow   stat
## Binarized data set
rev_b <- rev_c
rev_b$targV <- (rev_b$targV + 1) / 2
rev_b$maskV <- (rev_b$maskV + 1) / 2
rev_b$targEcc <- (rev_b$targEcc + 1) / 2
head(rev_b)
##      subj dom stairStart maskV targV targEcc   V2 rev_n rev_val targEcc_ targV_ maskV_
## 2161    1   1         -1     0   0.0     0.0  0.0     5    -0.9  central   stat   stat
## 2162    1   1         -1     0   0.5     0.0 -0.5     5    -0.9  central   slow   stat
## 2163    1   1         -1     0   1.0     0.0 -0.5     5    -0.7  central   fast   stat
## 2164    1   1         -1     0   0.0     0.5 -0.5     5    -0.9     near   stat   stat
## 2165    1   1         -1     0   1.0     0.5  0.0     5    -1.0     near   fast   stat
## 2166    1   1         -1     0   0.5     0.5 -0.5     5    -1.0     near   slow   stat
## A strictly linear model: will not capture nonlinearities in targEcc, targV & maskV
pvalfn(lmer(rev_val ~ dom + stairStart + rev_n + maskV * targV * targEcc + (1|subj), data=rev_c))
##                       estm    se       df    tVal  pVal star
## (Intercept)         -0.755 0.061   10.561 -12.359 0.000  ***
## dom                  0.040 0.057    8.000   0.700 0.504     
## stairStart           0.055 0.004 2681.001  12.436 0.000  ***
## rev_n                0.002 0.003 2681.001   0.769 0.442     
## maskV               -0.002 0.005 2681.001  -0.450 0.653     
## targV               -0.006 0.006 2681.001  -1.125 0.261     
## targEcc             -0.004 0.005 2681.001  -0.801 0.423     
## maskV:targV          0.031 0.007 2681.001   4.474 0.000  ***
## maskV:targEcc       -0.019 0.007 2681.001  -2.856 0.004   **
## targV:targEcc        0.003 0.007 2681.001   0.504 0.615     
## maskV:targV:targEcc -0.005 0.009 2681.001  -0.583 0.560
pvalfn(lmer(rev_val ~ dom + stairStart + rev_n + maskV_ * targV_ * targEcc_ + (1|subj), data=rev_c))
##                                      estm    se       df    tVal  pVal star
## (Intercept)                        -0.679 0.064   12.514 -10.654 0.000  ***
## dom                                 0.040 0.057    8.000   0.700 0.504     
## stairStart                          0.053 0.004 2662.001  13.053 0.000  ***
## rev_n                               0.002 0.003 2662.001   0.850 0.395     
## maskV_slow                          0.163 0.030 2662.001   5.527 0.000  ***
## maskV_fast                         -0.131 0.030 2662.001  -4.419 0.000  ***
## targV_slow                         -0.266 0.030 2662.001  -9.018 0.000  ***
## targV_fast                         -0.087 0.030 2662.001  -2.933 0.003   **
## targEcc_near                       -0.055 0.030 2662.001  -1.856 0.064    .
## targEcc_far                         0.004 0.030 2662.001   0.139 0.890     
## maskV_slow:targV_slow              -0.087 0.042 2662.001  -2.078 0.038    *
## maskV_fast:targV_slow               0.421 0.042 2662.001  10.076 0.000  ***
## maskV_slow:targV_fast               0.017 0.042 2662.001   0.405 0.686     
## maskV_fast:targV_fast               0.140 0.042 2662.001   3.343 0.001  ***
## maskV_slow:targEcc_near            -0.091 0.042 2662.001  -2.179 0.029    *
## maskV_fast:targEcc_near            -0.074 0.042 2662.001  -1.784 0.075    .
## maskV_slow:targEcc_far             -0.132 0.042 2662.001  -3.163 0.002   **
## maskV_fast:targEcc_far             -0.022 0.042 2662.001  -0.527 0.598     
## targV_slow:targEcc_near             0.089 0.042 2662.001   2.126 0.034    *
## targV_fast:targEcc_near             0.045 0.042 2662.001   1.078 0.281     
## targV_slow:targEcc_far              0.112 0.042 2662.001   2.672 0.008   **
## targV_fast:targEcc_far              0.028 0.042 2662.001   0.666 0.506     
## maskV_slow:targV_slow:targEcc_near  0.093 0.059 2662.001   1.573 0.116     
## maskV_fast:targV_slow:targEcc_near -0.087 0.059 2662.001  -1.476 0.140     
## maskV_slow:targV_fast:targEcc_near  0.032 0.059 2662.001   0.535 0.593     
## maskV_fast:targV_fast:targEcc_near -0.011 0.059 2662.001  -0.179 0.858     
## maskV_slow:targV_slow:targEcc_far   0.094 0.059 2662.001   1.597 0.110     
## maskV_fast:targV_slow:targEcc_far  -0.155 0.058 2662.001  -2.681 0.007   **
## maskV_slow:targV_fast:targEcc_far  -0.002 0.059 2662.001  -0.037 0.970     
## maskV_fast:targV_fast:targEcc_far  -0.070 0.063 2662.001  -1.121 0.262